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Aims. We investigate the spatial and temporal evolution of the heating of the corona of a cool star such as our Sun in a three- 
dimensional magneto-hydrodynamic (3D MHD) model. 

Methods. We solve the 3D MHD problem numerically in a box representing part of the (solar) corona. The energy balance includes 
Spitzer heat conduction along the magnetic field and optically thin radiative losses. The self-consistent heating mechanism is based 
on the braiding of magnetic field lines rooted in the convective photosphere. Magnetic stress induced by photospheric motions leads 
to currents in the atmosphere that heat the corona through Ohmic dissipation. 

Results. While the horizontally averaged quantities, such as heating rate, temperature, or density, are relatively constant in time, 
the simulated corona is highly variable and dynamic, on average reaching the temperatures and densities found in observations. The 
strongest heating per particle is found in the transition region from the chromosphere to the corona. The heating is concentrated in 
current sheets roughly aligned with the magnetic field and is transient in time and space. This supports the idea that numerous small 
heating events heat the corona, often referred to as nanoflares. 

Key words. Sun: corona — Stars: coronae — Magnetohydrodynamics (MHD) — Methods: numerical 
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1. Introduction 

The nature of the heating mechanism leading to a several million 
Kelvin hot outer atmosphere of cool stars still remains elusive. 
It is generally agreed that the mechanism heating the corona is 
related to the conversion from magnetic to thermal energy. One 
fundamental problem is that the actual dissipation of the (mag- 
netic) energy will appear on microscopic scales, while the ob- 
servable structures in the corona are macroscopic. A compre- 
hensive overview of coronal heating mecha nisms can be found 
in the classical conference proceedings of lUlmschneider et~aD 
(Il99ll) . For example, if we assume that a magnetic resistivity 
follows from classical transport theory, the dissipation should 
occur well below length scales of 1 m. Likewise, the ion (elec- 
tron) gyro-radii are close to m (cm) for typical coronal magnetic 
fields, so that the ultimate dissipation process has to operate on 
these small scales. 

In contrast to this, the actually observed coronal structures on 
the real Sun (and presumably on other stars) are on much larger 
scales. They range from several Mm above magnetic patches of 
the chromo spheric network to hundreds of Mm for large active 
region loops, i.e. are almost close to the solar radius. It is clear 
that the gap in length scale to the microscopic processes of a 
factor of some 10^ to 10^ cannot be bridged by a model that 
aims to describe the whole range of coronal structures. 

As a result, a large-scale model for actually observable 
solar structures has to employ a ( more or less sophis ticated) 
parametrization of the heating rate. iRosner et"ar 1 (Il978ll solved 
static one-dimensional (ID) models for coronal loops using dif- 
ferent parametrization, including the assumption of a (volumet- 
ric) heating rate that is constant in space and time, and de- 
rived their now classical scaling laws relating coronal temper- 
ature and pressure to heating rate and loop length. In a dy- 
namic ID loop model. iHansteenl (1 19931) assumed a heating rate 



that is transient in time and space by increasing the internal 
energy from one time step to the next at one grid cell in the 
numerical model, in order to study the response of the tran- 
siti on region to n anoflare-like heating as o riginally proposed 
by iParkerl (Il988b . iMiiller et al.l (l2003l l2004 show that a (vol- 
umetric) constant heating rate can lead to dynamic evolution 
and catastrophic cooling within a corona l loop if the h eating 
rate decreases exponentially with height. lAiouaz et~ar ] (I2005h 
investigated the consequences of diflTerent parametrization of 
the heating r ate on the outflow profile from a coronal funnel. 
[Patsourakos & KlimchuS (l2006l) modeled a loop composed of 
many individual (ID modeled) strands, each heated impulsively 
but constant in space. They find significant deviations from a 
Gaussian line profile of the emergent emission line, which still 
needs to be confirmed by observations. IWarren et al.l (1201 Ol) em- 
ploy multi- stranded loop models, too, by varying the energy in- 
put and the timing of the individual heating events in order to 
match observed spatial and temporal variations of the emission 
from coronal loops. In a more global approach, ISchrijver et al.l 
(120041) assume a parametrization as a function of the magnetic 
field and loop length for (approximate) static loop models. By 
comparing the appearance of the model coronae to the observed 
structure, they infer the details of this parametrization. 

This list of studies has to be incomplete, but all these (mostly 
ID) models have in common that they assume the parametriza- 
tion ad-hoc. By their very natur these models always suflfer from 
ignoring the spatial complexity and the changing 3D structure 
of the magnetic field in which the loop under investigation is 
embedded. Of course the major advantage of a ID model is the 
high spatial resolution that can be achieved and the possibility of 
including non-equilibrium ionization (e.g.. iBradshaw & MasonI 
2003ha). 

A 3D magneto hydrodynamic (MHD) model can account for 
the spatial complexity on the real Sun, and more important is 
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that it allows a more self-consistent tre atment of how the heat- 
ing rate is changing in space and time. iMok et af] (1201 Ol) use a 
parametrization depending on magnetic field strength. Through 
this they could compare the appearance of the loops in the 3D 
model to ID models. 

Th e 3D M HD coronal models by iGudiksen & Nordlundl 
(120021 l2005allbh include a heating based on field line braid- 
ing: foot point motions in the photosphere deform the mag- 
netic field, which induces currents that are subsequently dissi- 
pated. This Ohmic heating is intermittent in time and space as 
it depends on the evolution of the magnetic field which is self- 
consistently modeled. While the evolution of the magnetic field 
is treated self-consistently in the 3D MHD model, the Ohmic 
heating rate, rif should still be considered as a parametriza- 
tion. While the currents, j, are derived from the magnetic field 
in the numerical models, the magnetic resistivity, 77, used in the 
model is much greater than the valii e derived from transport the- 
ory (e.g. iBoyd & Sandersonll2003l) . This is due to limitations 
in the magnetic Reynolds number in numerical simulations. 
Therefore it is not clear whether 77 f represents the true heating 
rate. Considering results of the models of Gudiksen & Nordlund 
and the good match of these mo dels to solar spectroscopic obser- 
vations (^Peter et al l2QQ4ll2006h . it seems safe to at least consider 
rif as a (good) parametrization of the actual heating rate. 

This paper fo ll ows the m odel philosophy of 
IGudiksen & Nordlundl (l2002l l2005allbh and investigates the 
spatial and temporal variation of the heating rate in an active 
region. Special attention is paid to the energy distribution 
of individual energy releases, which can be considered as 
nanoflares. The outline of the paper is as follows. In Sect. |2] 
we describe our model and the set of MHD equations used. 
In Sect. [3] we analyze the data before we discuss the results in 
Sect. |4] and conclude in Sect.[5l 

2. The model corona 

The numerical model includes the solar atmosphere above a 
small active region in a 3D box of the size 50x50x30 Mm, which 
corresponds to roughly 0.05 solar radii. The dynamics and heat- 
ing of the corona stem from photospheric motions that braid the 
magnetic fields lines, often called flux braiding. Q This field line 
braiding essentially depends on the field geometry and the driv- 
ing boundary. Successively currents are induced that heat the at- 
mosphere by their dissipation. More precisely the photospheric 
motions, together with the magnetic field at high plasma beta, 
produce a Poynting flux into the upper atmosphere. The energy 
difl'erence between the non-force-free field and a potential field 
configuration, the free energy, is partly transferred to heat via 
dissipation of currents. The time scale of the conversion is given 
by the resistivity of the plasma. 

The idea of heat ing by Ohmic dissipation has already been 
proposed by iiParkery(ll983l) . He estimated the energy flux above 
an active region to be on the order of lOOWm"^ by comput- 
ing the magnetic stress introduced by the photospheric motions. 
The stress could be explained by a strain of magnetic lines of 
forces, which are then rapidly dissipated by reconnection. This 
paper follows that idea and investigates the rate and the spatial 
distribution of the dissipation. 

The temperature structure in the solar atmosphere is highly 
sensitive to th e radiative loss and the Spitzer heat conduction 
(ISpitzejll962l) . The latter is proportional to T^'^ and acts as a 

^ Strictly speaking, flux cannot be braided, thus this process should 
be called fleld line braiding. 
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Fig. 1. Initial vertical magnetic field component at lower bound- 
ary. 



thermostat for the corona. Energy input into the corona is con- 
ducted downwards along magnetic lines of force into regions 
of higher density. With increasing density at lower heights, the 
radiative loss becomes more efficient. If the heat input into the 
corona is increased, more energy is transferred into the chromo- 
sphere, and more energy has to be radiated away. Therefore the 
height of the transition region where the Spitzer heat conduc- 
tion brakes away moves down to heights of higher densities to 
compensate for the increased energy flux. As a result, the coro- 
nal density is higher when the atmosphere is in a pressure equi- 
librium. But the average coronal temperature change is small. 
Doubling the heating rate would only increase the temperature 
by a factor of 2^^^ due to the efficient heat conduction. In our 
model the energy input by the Poynting flux at photospheric 
level is dissipated into heat. Heat conduction transfers energy 
down to regions of high densities in the lower atmosphere where 
it is radiate d. We use the optically thin radiative losses given by 
ICook etal.1 ([19891 

The correct treatment of the energy equation is thus impor- 
tant for estimating the position of the transition region, as well 
as temperature and density of the corona. Direct c omparisons 
with observations with synthesized emission lines (iPeter et al.l 
i2004) are then possible. The intensity of the optically thin coro- 
nal emission lines is proportional to the density squared. Small 
changes in coronal densities dire ctly influ ence the coronal emis- 
sivities. In the model of Gudiks en & Nor dlund ( 2005a), the co- 
efficient of the Spitzer heat conduction is reduced by a factor of 
three. The heat conduction is the dominant process in the nu- 
merical scheme. Lowering the conduction results in larger time 
steps, thereby allowing for longer time series. Nevertheless, we 
use the coefficient as given in S pitzer & H arm ( 1953) and Spitzer 
(1962), i.e. three times more than in IGudiksen & Nordlund 
(i2005an . 

2.1. MHD equations 

The temporal evolution of the plasma and the magnetic field is 
governed by a set of partial diff'erential equations. These mag- 
neto hydrodynamic (MHD) equations are written in terms of the 
logarithmic temperature and logarithmic density, along with us- 
ing the vector potential for the induction equation. The former 
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is due to the huge variation from the photosphere to the corona, 
whereas the vector potential assures a solenoidal magnetic field 
at all times. The equations for the mass density p, the fluid ve- 
locity M, and the temperature T read as 

D\np 



Dt 



+ Vm = 0, 



(1) 



Du 1 / / \\ 

_ = -[-Vp+pg + jxB + 2vVo[pS_)) 



Dt 
DlnT 
Dt 



(2) 



+ (r - 1) V . w = Uo/ + 2pyS^ + + TVl . 

cvpT L — J 



(3) 

Because of gauge invariance, we can add the gradient of an ar- 
bitrary scalar field cp to the vector potential without changing the 
magnetic field 

5 - V X (A -h V0) . (4) 

We choose the resistive gauge cp = 77V • A. For constant 7/, the 
induction equation in terms of the vector potential reads as 



dA 

— = w X (V X A) 7/V^A . 
dt 



(5) 



The resistive gauge leads to a diff'usion of the vector potential 
proportional to the resistivity 7/, which is preferable in numerical 
simulations. Additionally the equation of state for an ideal gas 
correlates the temperature with the pressure, 

kB 



llMp 



-pT 



(6) 



We use the convective time derivative DIDt = d/dt -\- u • V to 
simplify the equations. The adiabatic constant y = Cp/cy is 5/2, 
the mean atomic weight /d is 0.667 and and ntp denote the 
Boltzmann constant and the proton mass. The gravitational ac- 
celeration is g = 274ms"^ and the viscous force is given by 
the gradient of the traceless rate of strain tensor S_, which only 
depends on derivatives of the velocity. 

The resistivity t] and dynamic viscosity v are constant and 
chosen so that the corresponding grid Reynolds numbers, i.e. R= 
Urmsdxv'^ and R= Wnnsdx7/"\ where dx is the grid spacing and 
Urms the rms velocity, in the model are close to one. For our grid 
resolution of several hundred of kilometers, we chose therefore 
rj = lO^^m^s and v = lO^^m^s. We apply a Newton cooling 
to the lowermost part of the model to adjust the temperature to 
follow a profile similar to the average model of IVemazza et al.l 
([T98lh . 

N=^(To-T), (7) 

^cool 

where T is the temperature, Tq = To(z) is the initial temperature 
profile, and the cooling time scale is given by 

Tcooi = Toexp(-z//z), (8) 

where z is the height in the box. The coefficients tq = 10"^ s and 
h = 40 km are chosen such that the influence on the temperature 
above 3 Mm is several orders of magnitude less than the other 
physical processes described in equation (|3]). Heat is transferred 
by anisotropic Spitzer heat conduction. The heat flux vector is 



(9) 



where ^0 = 10"^^ W (mK)"^ is the Spitzer value and b the unit 
vector of the magnetic field. For numerical stability we include 
mass diffusion on the lefthand side of equation ([T]) and isotropic 
heat flux proportional to |V In r|Vr into the heat flux vector in 
equation Q. 




40 60 80 
time t [min] 



Fig. 2. Averaged heating rate over time starting at t=0, e.i. initial 
conditions at t=0. 



2.2. Initial conditions 



iGudiksen & Nordlundl (l2005allbb used an observed magne- 
togram with an spatial extend of 250x250 Mm^. Because of 
the requested grid size, the magnetogram was scaled down by 
roughly a factor of 5 to fit into the computational domain. A res- 
olution of 400 km is required to resolve both the granular motion 
and the temperature gradient in the transition region. This down- 
scaling removes small-scale n iagnetic network patches. We u se 
the same magnetogram as in IGudiksen & Nordlundl (l2005allbh 
but to investigate the inffuence of the interaction between the 
active region and quiet Sun network fields small-scale features 
have to be introduced. The network ffux is taken from a second 
set of observations and is enhanced by a factor of five to in- 
crease the interaction with the active region. The resulting mag- 
netogram represents a small active region with a spatial extend 
of 50x50 Mm^ and is depicted in Fig. [T] To fill the computa- 
tional domain with magnetic ffux, we extrapolated a potential 
ffeld. The resulting corona was dominated by a large-scale loop 
connecting the main polarities of the active region. 

The initial plane-parallel temper ature and densi t y strat iffca- 
tion match the model atmosphere by lVemazza et al.l (Il98ll) . The 
atmosphere is at rest at the beginning. 

Once the simulation is started, after some 30 min solar time, 
the solution is independent of the initial setup; e.g., the typical 
granule life time is 5 min and the Alfven crossing time about one 
minute. Since we start with a potential ffeld, the averaged heating 
rate is zero at the beginning. After initialization time, the heating 
rate reaches values that are equal to the temporal average for the 
rest of the simulation time. Figure |2] shows the volume average 
of the squared current density beginning with the initial condi- 
tions. The dashed line in ffg.[2]marks the time after which data 
was taken for the analysis presented in this paper. This guaranties 
that data taken do not depend on the initial condition. The atmo- 
sphere becomes highly structured and dynamic with high veloc- 
ities, but the overall and averaged appearance remains more or 
less constant. 



2.3. Boundary conditions 

Photospheric motions at the lower boundary are the driving 
mechanism of the dynamics in the upper atmosphere. These mo- 
tions with power spectra similar to observed velocity spectra 
shuffle around the foot points of magnetic ffeld lines. This leads 
to non force-free magnetic ffelds with currents which heat the 
plasma by dissipation. The horizontal photospheric motions at 
the lower boundary are computed in a similar way as done by 
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iGudiksen & Nordlundl (I2002L 12005 allbh . These random-like hor- 
izontal motions are time-dependent and prescribe a velocity field 
similar to that found in the real solar photosphere from observa- 
tions. The vertical velocity at the bottom boundary, as well as all 
three components of the velocity vector at the top boundary, are 
zero at all times. The temperature T and density p at the bottom 
boundary are kept at their initial values. At the top boundary, the 
temperature and density have vanishing first derivatives in the 
vertical direction, but no specific values for T or p are imposed. 

The initial magnetic field configuration at the lower bound- 
ary evolves in time following the induction equation. Because 
photospheric random motions would corrode the active region, 
we update the magnetic field at the lower boundary by its initial 
value. The time scale involved is chosen so that the time series 
of the computed active region looks similar to an evolution of 
an active region on the Sun. At the top boundary we assume a 
potential field. 

In the horizontal directions the box is fully periodic. 
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Fig. 3. Horizontally averaged temperature (solid) and density 
(dashed) over height for one snapshot. Vertical dotted lines mark 
the magnetic transition region (at 4 Mm, cf . Sect. 13.1.11 ) and the 
maximum heating per particle (at 7 Mm, cf. Sect. 13.1.21 ). 



2.4. Numerical setup 

We use the Pencil Code (iBrandenburg & Doble j |2Q02|) to run 
our numerical model. It is a highly modular compressible MHD 
code tested on several astrophysical problems and can be used 
for massive parallel computing using Message Passing Interface 
(MPI). The numerical scheme is a finite diff'erence scheme com- 
prising a sixth-order spatial derivative and a third-order Runge- 
Kutta time- stepping scheme. The calculation is performed using 
128 grid points in each direction. This is the minimum amount 
of grid points needed to resolve the granular motions on an 
50x50 Mm^ area. The resulting grid spacing is on the order of 
390 km in the horizontal and 230 km in the vertical direction, 
respectively. 

Parameters are set to their values found in literature except 
the resistivity 77 and the dynamic viscosity v. Both are set to fulfill 
the requirement of a grid Reynolds number close to unity. The 
simulations runs for one solar hour before we start to collect 
data with a cadence of 30 seconds for another solar hour. This 
gives us 120 snapshots of all physical variables such as velocity, 
temperature, density, and magnetic field. 

The model is conducted on a cluster consisting of 32 Intel 
CoreDuo (TM) Xeon processors. The cluster is located at the 
Kiepenheuer-Institut for Solar physics, Freiburg, Germany. The 
typical time step of the simulation is about one millisecond solar 
time, which results in an overall computing time on the order of 
25.000 cpu hours. 

3. Results 

The model shows a highly dynamic and structured corona as 
a response to the photospheric driving imposed at the lower 
boundary. The total energy is balanced and the model corona 
reaches temperatures above one million Kelvin. Horizontally av- 
eraged temperature and density profiles (Fig. |3) show the diff'er- 
ent layers starting from the photosphere, the chromosphere, and 
the transition region followed b y the corona. The resu lting av- 
erage atmosphere is similar to a IVemazza et al.l (Il98ll) standard 
model. The density drops by several orders of magnitude before 
it has a large constant scale height above 8 Mm. The average 
coronal density is on the order of 10"^^kgm"^ or the particle 
number density lO^^m"^. 

The transition region in the averaged temperature profile 
spans a height of 2 Mm. Figure |4] shows the isosurface at lO^K 



above the vertical magnetic field component at the lower bound- 
ary. This highly corrugated surface represents the center of 
the transition region. The density along the isosurface is not 
constant, and therefore the emissivities of emission lines will 
vary. The transition region above the main polarities is at lower 
heights as the average transition region height. This is a re- 
sult of the anisotropic heat conduction. More heat is channeled 
along the magnetic field lines towards the main polarities than to 
the network flux because the corona is mostly connected to the 
strong magnetic flux concentrations. Thus the transition region 
migrates downwards to regions of higher density. The stronger 
heat income is then compensated for by increased radiative loss 
(cf. Sect.O. 

Figure |4] shows loop like structures with dense regions sep- 
arated by loops with lower density. The loops are filled and 
release their mass during the one-hour simulation. The pic- 
ture shows only one snapshot. A temporal analysis is done in 
Sect. [33] 

In Fig. [3] the magnetic transition region is marked at a height 
of 4 Mm. This height depends on the typical length scale of 
the magnetic features in the photosphere and is discussed in 
Sect. l3.1.T] The height of the maximum of the horizontally aver- 
aged heating rate per particle (cf . Sect. 13.1.2] ) is located at 7 Mm 
(Fig. [I. 

3.1. Average spatial distribution of the heating rate 

The average heating rate drops dramatically with height up to 
some 4 Mm and then falls off' more slowly but still exponentially 
(cf. Fig.O. Below 4 Mm most of the small-scale magnetic field 
lines emerging from the photosphere are closed back to the sur- 
face and thus causing the bend seen in Fig. [5] which is discussed 
in Sect. 13.1.11 In contrast, the average heating rate per particle 
peaks at some 7 Mm height (cf. Sect. 13.1.21) . Overall, the av- 
erage heating rate through Ohmic dissipation corresponds well 
with the average Poynting ffux into the upper atmosphere at all 
heights (cf. Sect. lTOl) . 

3.1 .1 . The magnetic transition region 

The vertical magnetic field at the lower boundary consists of two 
strong opposite polarities and small-scale network ffux patches. 
Field lines emerging from the active region reach high up into 
the atmosphere and build the coronal loop. Field lines starting at 
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Fig. 4. Visible impression of transition region height. Isosurface 
of the temperature at log r/[K]=5.0 and a vertical cut through the 
domain showing the density above the main magnetic polarities 
(gray-scale bottom picture). Color code of the isosurface and the 
plane indicates logarithmic densities. A dense coronal loop is 
visible in the vertical cut. The domain shown is 50x50x30 Mm. 
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Fig. 5. Horizontally averaged heating rate per unit volume for 
one snapshot. In the upper part above 4 Mm the average heating 
rate drops exponentially with an almost constant scale height of 
about 5 Mm. Vertical dotted lines as in Fig.O 



lower flux concentrations close back earlier and thus are shorter. 
The top panel of Fig.[6]shows a histogram of lengths of field lines 
traced from the lower boundary. The 256x256 starting points are 
equally distributed in the x-y-plane. The number of field lines 
decreases in roughly inverse proportion to length. But there is 
a change in the distribution function at roughly 15 Mm. Below, 
i.e. from to 15 Mm, the number of field lines per length inter- 
val decrease more slowly than exponentially. The total fraction 
of magnetic field lines in this interval is 75%. These field lines 
connect network patches and reach a height of roughly 4-5 Mm 
when semi-circular loops are assumed. For lengths above 15 Mm 
the number of field lines decreases roughly exponentially up to 
some 40 Mm. For even longer field lines we find a clustering 
with the longest field lines reaching almost the top of our do- 
main. These field lines connect the main polarities and are less 
than l%in total. 

A distinct change in the distribution of field line lengths is 
visible in Fig. [6] at a length of about 15 Mm. This field line 
length corresponds to an apex height of about 4 Mm (assum- 
ing a semi-circular loop). This indicates that at heights in the 
computational domain below roughly 4 Mm the magnetic field 
topology is dominated by short loops connecting small-scale fea- 
tures. Above this height the volume is dominated by the longer 
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Fig. 6. Top panel: Histogram shows the distribution of field line 
length in the domain respecting the periodic boundaries. Equally 
distributed at z=0 have been selected for the tracing algorithm 
256^ starting points. Vertical dotted line depicts the local maxi- 
mum at 15 Mm. 



field lines connecting the two main magnetic patches of the ac- 
tive region. We therefore denote this height range at about 4 Mm 
where the magnetic topology changes from small-scale to large- 
scale as the magnetic transition region. This can also be seen as 
a kink in the distribution of currents (or more exactly the heat- 
ing rate o cp-) with height in Fig. [5] (Sect. l4.2l) . This definition is 
similar to I Jendersie & Peterl (120061) . 



3.1 .2. Foot-point dominated heating 

We analyzed the heating rate for each snapshot of our model by 
investigating the Ohmic heating rate (cf . Eq. ©) ^yUoj^ at all grid 
points. It corresponds to the volumetric heating rate and Fig. [5] 
shows its horizontal averages for one time step. The volumetric 
heating rate decreases over more than ten orders of magnitude 
from the photosphere to the corona. 

From the photosphere to the upper chromosphere at 4 Mm, 
the heating rate decreases about eight orders of magnitude. Thus 
99% of the energy is deposited in the chromosphere, and the 
chromosphere acts as an energy buff'er. 

Above 4 Mm (lefthand side in Fig. [5]) the heating rate drops 
exponentially with a scale heig ht being on average around 5 Mm. 
iGudiksen & Nordlundl (l2005bl) found an average scale height of 
5 Mm in their model, too. For both the transition region and 
the corona, the heating scale height is constant. The volumet- 
ric heating rate at the top of the chromosphere is approximately 
10-5 Wm-^ 

The density scale height (cf. Fig. [3]) at roughly lO'^K (below 
5 Mm) is about 0.3 Mm. Due the rapid temperature increase up 
to 10^ K the density scale height becomes 13 Mm above 7 Mm. 
Between the heights of 5 Mm and 7 Mm, the volumetric heating 
rate drops exponentially with a scale height between the density 
scales at these levels. This leads to a maximum specific heating 
rate per particle as illustrated in Fig. |7] The specific heating rate, 
i.e. per particle, increases starting at the lower chromosphere 
with a scale height of 0.5 Mm, whereas the specific heating rate 
drops exponentially with a scale height of 6 Mm above 7 Mm 
height. Thus, even though most of the energy is deposited in the 
chromosphere, the available energy per particle is higher in the 
corona. 

The volumetric heating rate (cf. Fig.O, as well as the heating 
rate per particle (cf. Fig. IT}, shows that the heating of the coronal 
loops is foot-point dominated, and yet there is heating also in the 
upper part of the corona. 
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vertical z [Mm] 

Fig. 7. Heating rate per unit mass. Scatter plot shows the distri- 
bution of logarithmic specific heating rate. Overplotted is log- 
arithm of the horizontally averaged specific heating rate, i.e., 
heating per particle, for one snapshot. Vertical dotted lines as 
in Fig.O 
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vertical z [Mm] 

Fig. 8. Horizontal, averaged, upwardly directed heat flux (solid) 
derived from volumetric heating rate. Horizontally averaged ab- 
solute Poynting flux is overplotted as a dashed line. Both are 
derived from one snapshot in time. Vertical dotted lines as in 
Fig. [3] 



even valid for all heights. The energy released in a sub- volume 
corresponds to the incoming Poynting flux. 



3.1 .3. Energy flux into the atmosphere 

The volumetric heating rate can be converted into an energy flux 
by assuming that all energy comes from the lower boundary. 
This is the case for the Ohmic heating as the driving occurs at 
the photosphere and the magnetic field dominates the upper at- 
mosphere. The energy flux density 2 is a function of height and 
gives a measure for the energy per unit time that has to pass this 
height through a unit plane. The value of Q can be derived from 
the volumetric heating rate by integrating from z to infinity; i.e., 
the heating above z is powered by the energy flux through the 
height z, 

foo 
wd(x,y.zf^z . (10) 

The horizontal average of the energy flux density is depicted 
in Fig. [8l The energy flux into the upper atmosphere is ranges 
from 10^ to 10^ Wm"^ and decreases rapidly up to 4 Mm in 
height, the magnetic transition region. At the bottom of the tran- 
sition region the energy flux density is a few times 100 Wm"^. 
This is consistent with the typical observation-based estimations 
for the energy flux need ed to heat and sustain the corona (e.g. 
IWithbroe & Noveslll977h . The energy flux decreases further to a 
few times 10 Wm"^ in the upper part of the corona. 

The Poynting vector describes the direction and the strength 
of the energy flux in an electro-magnetic field. Using the induc- 
tion equation © and Ohm's law, we can write the Poynting vec- 
tor S independent of the electric field as 

S= —(woj-uxB)xB. (11) 

Figure[8]depicts the horizontally averaged vertical component of 
the Poynting vector. As the Ohmic heat is energy converted from 
the magnetic field, the Poynting flux roughly follows the energy 
flux derived from the heating rate. Because the Poynting flux 
depen ds on t he plasma velocities it shows more temporal vari- 
ation. iGalsga ard & Nordlund (1996) found that the dissipation 
rate scales linear with the Poynting flux at the lower boundary. 
Averaging over a short time period we find that this relation is 



3.2. Current sheets and nanoflares 

The horizontally averaged heating rate in the simulation is al- 
most constant in time. Investigating smaller volumes, i.e. at high 
spatial resolution, the heating rate reveals a highly dynamic na- 
ture with a broad range of time scales. Figure |9] depicts heating 
rate in a slab placed above the main polarities of the magne- 
togram. The heating rate is normalized by its horizontal average 
in order to remove the steep vertical gradient (cff' Fig. [5]). 

We used a narrow 3D volume, i.e. a slab, rather than a 2D 
plane and integrate over one direction. A vertical cut through 
the domain would only show the intersection of the magnetic 
field lines with the plane but not a part of the loop. Samples of 
field lines that fit completely into the slab are shown in Fig. \TU\ 
The field lines seem to intersect, which is only an eff'ect of the 
projection onto the x-z plane. 

In the coronal part the Ohmic heating rate is organized in cur- 
rent sheets oriented parallel to magnetic field lines. These struc- 
tures exist at the resolution limit and are several grid cells wide. 
The alignment of the heating structures to the field lines is illus- 
trated in Fig. [TOl Below 4 Mm, the magnetic transition region, 
the normalized heating rate shows small structures that indicate 
the short field lines (cf. Sect. 13. Lit . 

The normalized Ohmic heating rate in Fig. [9] only illustrates 
the spatial distribution of the heating but does not give any mea- 
sures. The specific heating rate per particle in the same slab is 
depicted in Fig.[TT] It also shows the concentration of the heat- 
ing in structures along the magnetic field lines. Furthermore, the 
foot-point dominated heating is illustrated. As shown in Fig. [7] 
the specific heating rate peaks at around 7 Mm. 

Current sheets along magnetic loops have also been investi- 
gate by e.g. iRappazzo et al.] (l2007l) in high resolution loop sim- 
ulations investigating the Parker field line tangling that leads to 
Ohmic dissipation. Their model compares well to our setup in- 
cluding constant resistivity, whereas they used a so-called hyper- 
resistivity. But in contrast our model extends over larger volume 
and comprises a wide range of magnetic structures in a more 
realistic geometry. 
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Fig. 9. Integrated heating rate in a slab of the 
size Ax=50Mm, Aj=7Mm and Az=30Mm. 
Color coded is the heating rate averaged along 
the y direction divided by its horizontal mean. 
The figure is periodic in the horizontal direc- 
tion. Horizontal dotted lines as in Fig. [3] 
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Fig. 10. Integrated heating rate as in Fig. |9] 
but overplotted are the magnetic field lines in 
the slab projected onto the x-z plane. Starting 
points for the tracing algorithm are equally dis- 
tributed in the 3D slab. 
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Fig. 11. Logarithmic heating rate per particle in 
the same slab as in Fig. |9] averaged in the y- 
direction. Horizontal dotted lines as in Fig. |3] 



3.3. Heating dynamics and response along field lines 

In the previous section we showed that the volumetric heating 
rate is organized in current sheets. A closer look at the cur- 
rent sheets reveals that these current sheet are discontinuous. 
Furthermore, the discontinuities change in time. Therefore we 
analyzed the time dependence of the heating rate along magnetic 
field lines. Figure [12] illustrates the selected six diflTerent field 
lines that are traced in the complete 3D box respecting the peri- 
odicity of the domain. These field lines are selected to represent 
diflTerent heights and diflferent connectivities between the main 
polarities, as well as into network patches. A list of properties of 
the selected field lines is given in Table [T] In the subsequent dis- 



cussion of the physical properties, we use the word loop instead 
of the mathematical construct of a field line. Loop # 2 reaches 
the bottom of the transition region, whereas loops # 1 and 3 ex- 
tends above the transition into the base of the corona. Loops # 4 
to 6 extend into the corona. In comparison to half circles, the 
apex height of loops # 1 to 3 is smaller than half their foot-point 
distance, implying that these loops appear somewhat flattened. 
Loops #4 to 6 are stretched out into the corona with heights ex- 
ceeding half the foot-point distance. 

For the following investigation we did not follow the field 
lines in time, instead study the plasma properties in time on 
the trajectory defined by the field line at the beginning of the 
time interval under study. This is justified by the low (e.i. some 
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Table 1. Properties of loops shown in Fig.[T2l 



No 


apex height 


length 


foot point distance 




[Mm] 


[Mm] 


[Mm] 


1 


7.1 


28.9 


22.8 


2 


5.0 


26.5 


22.9 


3 


8.9 


33.8 


25.2 


4 


24.7 


76.7 


25.0 


5 


23.7 


76.5 


26.7 


6 


29.0 


87.5 


48.9 




Fig. 12. 3D representation of six magnetic field lines used for 
further analysis. White lines indicate the boundaries for the do- 
main periodic in the horizontal directions. The box size is there- 
fore 100x1 00x30 IMm. 



km/s) velocities perpendicular to the magnetic field and the cor- 
responding displacement over 10 minutes is close to the grid 
spacing. 

3.3.1 . Specific heating rate per particle 

Figure [13] displays the specific heating rate per particle for one 
solar hour along the magnetic field lines. Diff'erent types of heat- 
ing events can be found. On the left side of panels [5] and [6], the 
heating is continuous over more than 20 minutes. On the other 
hand, loop # 3 shows at the top (middle of panel [3]) short-lived 
heating events. Since the numerical time step is only a few mil- 
liseconds these events are resolved well in the simulation. 

Loops # 1 and 2 connecting the main polarities with the net- 
work field are mainly foot point heated at photospheric levels; 
however, the heating of the foot points in the network (cf . right- 
hand side of panels [1] and [2]) is stronger than the heating at 
the negative main polarity of active region. As the velocities are 
quenched at strong magnetic fields, the shear of the foot points 
becomes larger for the network flux patches. 

Loop # 3 undergoes two main heating events at its top. The 
first event (from 15min to 30min) consists of several short 
small-scale events. The energy content of these small events 
can be computed and is similar to nanoflares, i.e. some 10^^ J, 
so that the first event compares to nanoflare heating or heating 
by nanoflare storms. The second (starting at 30min) event starts 
with nanoflare heating before the heating gets stronger. The en- 
tire top of loop # 3 is heated in the last couple of minutes. 

Loop # 4 is heated mainly at the height of the transition re- 
gion where the heating events have durations between 10 and 
20min, and they occur on both sides of the loop. Loops #5 and 
6 undergo a long heating event at their foot points (on the left- 
hand side of panels [5] and [6]). 

As already shown in Figs. [5] and [7] loops are predominantly 
heated at their foot points. However, exceptions, e.g. loops # 1 or 
3, can be found. 



3.3.2. Temperature variation 

The temporal evolution of the temperature is depicted in Fig.[T4l 
Loop # 1 partly increases in temperature by an order of a mag- 
nitude at the beginning of our time series. At this time also the 
specific heating rate shows a peak. 

Loop #2 is the coolest and lowest loop. Its maximum tem- 
perature at the top decreases with time, nevertheless, heating 
events seem to occur that are not clearly visible in panel [2] in 
Fig. \T3\ Viscous heating has to be taken into account to explain 
the small temperature variations of this loop. Because this loop 
is entirely embedded in the chromosphere, the temperature (in- 
ternal energy) is al so slightly influenced by the Newton cooling 
term(cf. SectlTTI). 

Loop # 3 shows a nice correlation between temperature vari- 
ation and the specific heating rate, but in comparison the temper- 
ature along the field line is less structured. Due to the high eflS- 
ciency of the anisotropic Spitzer heat conduction, the internal en- 
ergy is distributed along the field line on short time scales. Small 
temperature variations due to the nanoflare heating are smeared 
out rapidly. 

Loops # 4 to 6 also illustrate a nice correlation between the 
specific heating rate and the temperature. At times when the foot 
points are not heated the loop cools down. When the heating sets 
in not only small regions but also the entire loop gets hotter as 
a result of the efficient Spitzer heat conduction. For short time 
periods the loops are almost isothermal. 

3.3.3. Vertical velocities 

Figure [15] shows the temporal evolution of the vertical compo- 
nent of the plasma velocity parallel to the magnetic field lines. 
The velocity is divided by the local sound speed to resolve the 
large difference between the photospheric values of a few km s"^ 
and the coronal ffows up to some 100 km s"^ In general, the ve- 
locities in the box are subsonic. 

By normalizing one can qualitatively distinguish between the 
types of motions that are a response to the heating along the field 
line. The signature of the granular motions in the lower bound- 
ary is seen as periodic changes of up and down ffows at the foot 
points of the ffeld lines. These motions have a period of roughly 
5 min, which corresponds to the lifetime of the granulation driv- 
ing the magnetic ffeld in the photosphere. 

Loop # 1 is not heated much, and it cools down for the given 
time series. One effect is that the loop drains. The apex seems to 
move upwards, which can be understood as a motion of the ffeld 
line itself. As the motion is much slower than one kms"\ ffeld 
line tracing can be neglected, as argued above. 

Loop # 2 shows a mix of up ff ow and down ff ow events along 
the loop. A clear relation cannot be found when this is compared 
to the speciffc heating rates. The chaotic like heating events re- 
sult in quite irregular ffow patterns. 

Loops # 3 to 6 again show a nice correlation between the spe- 
ciffc heating rate (Fig. and the vertical velocities (Fig. [T5]) . 
At places of strong heating the plasma evaporates and is fflled 
into the loop. The velocities are upwardly directed (blue color) 
and ranges from 50kms"^ to over 100 km s"^. At times of no 
heating, the coronal plasma cools down by anisotropic heat con- 
duction and radiative losses and starts to drain (red color) out of 
the loop. 

Loop #3 reveals a strong downffow on the lefthand side. 
There the loop is much denser than the opposite foot point. Thus 
the speciffc heating per particle is less and the radiative cooling 
is more efficient. The loop cools down quickly, and the plasma 
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normalized loop length normalized loop length normalized loop length 



Fig. 13. Heating rate per particle for one solar hour along magnetic field lines depicted in Fig. [121 Loop lengths are given in Table [T] 
The color table varies for each field line and is given on the righthand side of each panel. Loops # 1 to 3 are short, reaching any 
height below 10 Mm, while loops #4 to 6 represent coronal loops (see Fig.[T2l). 




normalized loop length normalized loop length normalized loop length 



Fig. 14. Same as in Fig. \T3\ but showing logarithmic temperatures. The color table varies for each field line and is given on the 
righthand side of each panel. 

has to drain out of the loop. After some 45 min, the heating be- right. These siphon flows are a result of the asymmetric heating 
comes strong enough so that plasma evaporation can take place of the foot points. First the foot point on the righthand side is 
again. heated, then the one opposite it. After some 45 min, the loop is 

above one million degrees on both sides due to the heat conduc- 

Loop # 6 shows a siphon flow first from the right to the left 
side, and after some 15 min the siphon goes from the left to the 
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tion. The draining stops and at both sides the loop is filled with 
plasma. 

Even though loops #4 to 6 are foot-point heated either on 
one or the other, a siphon flow is suppressed as soon as the loop 
approaches isothermal. The highly efficient Spitzer heat con- 
duction leads to the almost constant temperatures for which the 
loops either are filled by plasma or they drain. 



pation process. As long as the Poynting flux into the domain 
is constant, th e dissipat ion ra te is fixed on larger time scales. 
iGalsgaard & N ordlund (1 1996b find a scaling law between the 
dissipation rate and the Poynting flux that only depends on the 
granular motions and the magnetic field in the photosphere. This 
justifies our choice of // to have a magnetic Reynolds number of 
order unity. 



4. Discussion 

Our approach of a 3D MHD numerical model successfully pro- 
duces a one million degree hot corona a bove a small active 
region. The heating mechanism following iParked (Il978l) self- 
consistently produces a heating rate strong enough to balance 
Spitzer heat conduction and radiative loss. The empirically based 
photospheric motions create the dynamics of the system. The 
model corona is in a quasi- stationary state where energy in- 
put by the Poynting flux is balanced by radiative loss. Other 
3D models also find a hot corona b ut use additio nal empiri- 
cally based heating mech anisms (e.^.|Abbettll2007l). or else the 
spatial extend is less (e.g. iMartmez-Sykora et al.ll2009l) than in 
th e model presented in this paper. In contrast to the models 
of [Gudiksen & Nordlund' (2002, 2005 a'D), which are similar to 
ours, we also included the magnetic flux of the network patches. 



4.1. Notes on the heating rate 

The heating rates found in our model are sufl&cient to sustain 
a hot corona. Implemented Ohmic heating is only an approxi- 
mation because the magnetic resistivity used in the simulation 
is much greater than one would expect for the coronal plasma 
following classical transport theory (e.g Spitzer 1962). The re- 
sistivity 77 describes the process of conversion of magnetic to 
thermal energy, and the value used in our model is increased by 
about eight orders of magnitude, in order to have magnetic grid 
Reynolds numbers of order unity. So that the currents are actu- 
ally dissipated. 

The magnetic Reynolds number describes the proportion of 
the advective to the resistive process for a certain length scale or 
the proportion of the resistive time scale to the advective time 
scale 

Urms^ tresis (12) 



Rm = 



^advec 



where Urms is the rms velocity and / a given length scale. If we 
assume coronal values for the resistivity (e.g. lSpitzerlll962h . then 
we have to look at cm length scales to obtain a Reynolds number 
of 1 where the currents are dissipated. These length scales are 
not resolved, but we must use the grid spacing to compute the 
magnetic grid Reynolds number. Using / =400 km results in a 
resistive time scale several orders of magnitude larger than the 
advective time scale. Thus we would not dissipate any currents in 
the time of our simulation. This contracts observations, beacause 
the life time of active regions or eruptive events such as flares 
illustrates that the resistive process is much faster. To circumvent 
this dilemma we chose the resistivity t] such that the magnetic 
grid Reynolds number is on the order of unity and therefore the 
time scales are comparable. Naturally we find a wide spread of 
Re ynolds numbers in our domain. 

[Gals gaard & Nordlundl (1 19961) also investigated the depen- 
dency of the resistivity on the numerical resolution. For the 
first event in their model they found a logarithmic scaling, but 
for the rest the dissipation seems to be independent of reso- 
lution. Lowering the resistivity would just postpone the dissi- 



4.2. The magnetic transition region 

At a height of about 4 Mm the magnetic structure of the atmo- 
sphere is abruptly changing (Sect. 13.1.11 ). This manifests itself 
most clearly in the kink of the average heating rate and the en- 
ergy flux as a function of height (Figs. [5] and (8]). This is related 
to the height where most of the small-scale (granulation scale) 
magnetic flux closes, defining a magnetic transition region. Only 
above this height does the upper atmosphere reach a magnetic 
state that represents the magnetic structure of the corona. 

This magnetic transi tion region is loc ated abo ve the heigh t 
of the classical canopy (lGiovanellilll980 : Solank i et al"]ll992h . 
which is found below 1 Mm. The latter is set by the rapid ex- 
pansion of the magnetic field with height, because the gas pres- 
sure drops exponentially and a horizontal equ ilibrium of gas and 
magnetic pressure has to be achieved (e.g. 'Solanki & Steinerl 
0,9 90) . Basically the classical canopy is a magneto-hydrostatic 
eff'ect. 

In contrast, the height of the magnetic transition region can 
be understood by the potential field extrapolation of the distribu- 
tion of magnetic flux at the surface. In a numerical experiment 
iJendersie & Petej (l2006l) show that, in quiet Sun regions above 
the network, the small-scale magnetic concentrations would pro- 
duce numerous short loops reaching up to about 4 Mm (cf . their 
Fig. 3). The small-scale structures push up the expanding field 
lines from the larger magnetic patches, so in contrast to the 
classical canopy, where the (predominantly horizontal) mag- 
netic field is overlying a field-free region (e.g. lSteiner & MurdinI 
[2000) a larger- scale magnetic field is found above a volume with 
small-scale closed magnetic field lines, in the case of the mag- 
netic transition region . 

Our current model carries this concept further by giving 
up the assumption that the magnetic field is potential in na- 
ture (i.e. current- free). The resulting currents found in the 3D 
MHD model are very strong in the photosphere and chromo- 
sphere because of the shear applied by the horizontal foot-point 
motions of the granular convection. The heating rate rjf drops 
roughly exponentially with a very small scale height of less than 
0.5 Mm. Above the magnetic transition region, where the larger 
scale magnetic structures dominate, the heating rate still drops 
exponentially (on average), but now much more slowly with a 
ten times greater scale height of about 5 Mm. 



4.3. The chromosphere as an energy filter 

While the heating is highly intermittent in time and space, 
the horizontally averaged heating rate is almost constant in 
time. At the magnetic transition region, i.e. at the top of the 
chromosphere, the volumetric heating rate drops to just be- 
low 10""^ Wm"^. There the energy flux into the corona heat- 
ing up the upper atmosphere is a bit more than 100 Wm"^. 
This is consistent with typical es timates of the energy dem ands 
derived from observations (e.g. IWithbroe & Noyesfll977l) . es- 
pecially when considering that we describe a small active 
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Fig. 15. Same as in Fig. [131 but showing vertical component of flow velocity parallel to magnetic field lines. The velocity is divided 
by the local sound speed. The color table varies for each field line and is given on the right hand side of each panel. Red colors 
indicate down flows and blue colors up flows accordingly. In general, flows are subsonic. 



region in our model. Our results are also consistent with 
iGudiksen & Nordlundl(l2002L 12005 allbh . 

It is most remarkable that in our 3D MHD model the up- 
wards directed energy flux heating the corona is dropping by 
about six orders of magnitude (Fig. [8]) from the surface to the 
base of the corona, i.e. the magnetic transition region. This is 
consistent with the real Sun, because in the model we have 
roughly the same Poynting flux at the surface as found on the 
Sun. This implies that only the 10"^th part of the energy flux at 
the surface makes it into the corona. Or in other words: about 
99.9999% of the energy is dissipated in the lower atmosphere 
below the magnetic transition region! It has to be stressed that 
no fine-tuning was applied to the model to get the correct energy 
flux into the corona, i.e. to dissipate just 99.9999% (and not all) 
of the energy in the lower atmosphere. 

Thus the lower atmosphere below the magnetic transition 
acts as an efl&cient filter for the energy transported upwards to 
heat the corona. It would be very interesting to investigate how 
the efficiency of this filter changes with the parameters of the 
3D MHD model and the boundary conditions; i.e., what can be 
expected for other types of stars than our Sun. This could af- 
fect how we understand why the X-ray l uminosity in more ac - 
tive stars is strongly enhanced (up to lO'^ HPizzolato et al.ll2003l) , 
while the filling factor of the coronal plasma can be increased by 
a factor of only about 100 (assuming a filling factor for the Sun 
of about 1%). Further studies will have to elucidate this problem. 

4.4. Individual strands and coronal loops 

While in a ID numerical model for a coronal loop one can aff'ord 
a high spatial resolution to resolve strong gradients and shocks 
along the loops or include ionization processes, a 3D MHD 
model provides the possibility of investigating how neighbor- 
ing structures interact and getting a better (more self-consistent) 



description of the heating rate. Therefore it is interesting to com- 
pare our 3D model to state-of-the-art ID loop models, as well as 
to multi- stranded loop models. 

4.4.1 . Thermal and dynamic response of individual loops 

If we investigate individual field lines in our 3D computational 
box, we can follow the temporal evolution of the physical pa- 
rameters in the same manner as in ID models for coronal loops 
(Sect. [33l Figs. [13] to [15]). This shows a thermal and dynamic 
evolution of the loops, i.e. variations in temperature and veloc- 
ity, on time scales from well below minutes to one hour. 

While the time scale of the driver in the photosphere is 
only several minutes (granulation), the response of the corona 
in terms of the heating rate shows much faster variations. This 
is a consequence of the non linearity of the physical process of 
field line braiding and energy conversion. Therefore the rapid 
variations in temperature and velocity seen in the coronal parts 
of individual loops are not a signature of the photospheric driv- 
ing, but the response of the plasma to the small heating events 
in rapid succession. In that sense the dynamics of the corona is 
decoupled from the photosphere. (In the lowermost parts of the 
loops a clear signal of the photospheric driver can be seen, of 
course, cf. Fig. [15]). 

The variations on longer time scales are (partly) due to the 
long-term evolution of the larger magnetic field patches in the 
photosphere. For example, over a half hour the granular motions 
can create (or destroy) by chance larger patches of stronger mag- 
netic field, which in turn result in an increased heating of the 
overlying corona (e.g. loops #5 and #6 in Fig. [13]). In general the 
evolution of the temperature of a loop on these longer time scales 
roughly follows the heat input from below, but with some modi- 
fication because of the long cooling time of the coronal plasma. 
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In general the dynamic evolution as found along field lines 
in our 3D models compares well with state-of-the-art ID loop 
models. For example, we find dynamic processes such as cool- 
ing, draining, or siphon flows of the same o rder of magni t ude as 
comp arable gradients as described, e.g., in iMiiller et al.l (l2003l 
|20Q4|) . While they assumed an exponential drop of the heating 
rate in their loops with a comparable scale height, they kept the 
heating constant in time. In the light of our finding that the de- 
tails of the photospheric driver (on short time scales) are not 
important for the coronal dynamics, this can be considered as 
a minor diff'erence. Therefore we can conclude that our results 
for individual loops forming in the 3D MHD model compare 
well to high-resolution ID loop models. Of course, the compu- 
tational eff'ort prevents us from including, e.g., ionization pro- 
cesses, but we can properly resolve the thermal and the flow 
structure within the loops. This allows us to also draw conclu- 
sions on multi- stranded loop structures. 

4.4.2. Comparison to multi-stranded loop models 



In the ir multi- stranded loop models iPatsourakos & KlimchukI 
(l2QQ6h assume that a coronal loop consists of individual strands 
that are treated as independent ID models. This is justified by the 
low heat conduction across the magnetic field, which basically 
isolates the individual structures. Their strands are heated by in- 
creasing the heating rate for some time (and by this mimicking 
nanoflares). In our 3D model a large loop, such as seen in the 
vertical cut of Fig.lH contains numerous loop-like current sheets 
that are parallel to the magnetic field, as depicted in Figs. [9] and 
[TOl These current sheets can be very thin, down to the size of the 
grid cells of the computation. Thus one could consider a bun- 
dle of field lines with the parallel current sheets (Fig. [TOll as the 
many strands making up a loop in the multi- strand model. 

However, there is a major diff'erence to the multi- 
strand model. While the cross-field conduction is very low, 
i.e. the strands are thermally isolated, the 3D model in- 
corporates interaction of the strands: the braiding of the 
magnetic field lines causes the loop- shaped current sheets 
parallel to the magnetic field. This results in a self- 
consistent impulsive heat i ng of the individual strands, whereas 
IPatsourakos & KlimchukI j2QQ6l) release the nanofl are energy 
ad hoc. Recently Lopez Fuentes & KlimchukI (|201Q|) have em- 
ployed a multi-strand model that mimics this process. They as- 
sume that the strands are displaced by photospheric motions, 
which set s the time for the energy release for a given strand ac- 
cording to lParkerl(ll988l) . but there is still no interaction between 
the stands as these are desc ribed by ID models. 

Another difl'erence to IPatsourakos & KlimchukI (l2006l) is 
that they assume the heating is distributed uniformly along the 
loop, while our model gives a heati ng rate that drops expo- 
nentia lly with height. However, as the IPatsourakos & KlimchukI 
(l2006l) model loops are much larger than those fitting in our com- 
putational box, it remains to be seen if this diff'erence is signifi- 
cant or if this is more a property of smaller than of larger ones. 

More analysis of the distribution of the heating rate in space 
and time as found in our 3D MHD model will have to show 
how this compares in detail with the assumptions of the multi- 
stranded loop model, and might provide important constraints 
for the further development of the multi- stranded loop models. 

5. Conclusion 

Our 3D numerical model of the solar corona successfully de- 
scribes how to sustain a hot corona by the heating mechanism 



based on Ohmic dissipation. The average heating rate and the 
derived energy flux compare w ell to the observational require- 
ments (e.g. Withbroe & Noyes 19771) . Furthermore, the model 
resolves the intermittent and transient character of the heating 
on a wide range of energy scales. 

The heating rate per particle (or mass) is found to be 
strongest in the transition region from the chromosphere to the 
corona. This is because of the scale height of the (on average) 
exponentially dropping volumetric heating rate is between the 
pressure scale heights of the chromosphere and of the corona. 

The heating is concentrated in current sheets that are roughly 
aligned with the magnetic field lines and which are highly inter- 
mittent in time and space. In general this supports the idea that 
the corona is heated by a large nu mber of smal l energy deposi- 
tions, often named nanoflares (e.g. |Parke3ll988h . 

The dynamics within a corona loop (or a strand thereof), i.e. 
single field lines in the complex 3D magnetic field, follow the 
transient heating events. Loops with siphon flows, as well as ir- 
regular flow patterns, can be found in the 3D MHD model. The 
time-dependent heating rate along individual field lines as found 
in our 3D MHD models may be used as input to higher resolu- 
tion ID loop models including, e.g., non-equilibrium ionization. 

The results show that 3D numerical box models of the corona 
are a useful tool to investigate the nature of coronal heating, in 
particular the distribution of the heating rate in time and space. 
The transient heating can be investigated down to energies well 
below energetic events currently observed on the Sun. 

Acknowledgements. We would like to thank Wolfgang Dobler for introducing us 
to the Pencil Code. Special thanks are due to Boris Gudiksen for the discussions 
when we started this project and for providing the code for the photospheric 
driver that we used in earlier experiments. This work was supported through 
grants by the Deutsche Forschungsgemeinschaft (DFG). 



References 

Abbett, W. P. 2007, ApJ, 665, 1469 

Aiouaz, T., Peter, H., & Keppens, R. 2005, A&A, 442, L35 

Boyd, T. J. M. & Sanderson, J. J. 2003, The Physics of Plasmas (Cambridge 

University Press) 
Bradshaw, S. & Mason, H. 2003a, A&A, 407, 1127 
Bradshaw, S. & Mason, H. 2003b, A&A, 401, 699 

Brandenburg, A. & Dobler, W. 2002, Computer Physics Communications, 147, 
471 

Cook, J. W., Cheng, C.-C, Jacobs, V. L., & Antiochos, S. K. 1989, ApJ, 338, 
1176 

Galsgaard, K. & Nordlund, A. 1996, J. Geophys. Res., 101, 13445 
Giovanelli, R. G. 1980, Sol. Phys., 68, 49 
Gudiksen, B. & Nordlund, A. 2002, ApJ, 572, LI 13 
Gudiksen, B. & Nordlund, A. 2005a, ApJ, 618, 1020 
Gudiksen, B. & Nordlund, A. 2005b, ApJ, 618, 1031 
Hansteen, V. H. 1993, ApJ, 402, 741 
Jendersie, S. & Peter, H. 2006, A&A, 460, 901 
Lopez Fuentes, M. C. & Klimchuk, J. A. 2010, ApJ, 719, 591 
Martmez-Sykora, J., Hansteen, V., DePontieu, B., & Carlsson, M. 2009, ApJ, 
701, 1569 

Mok, Y., Lionello, R., Mikic, Z., & Linker, J. 2010, in American Astronomical 
Society Meeting Abstracts, Vol. 216, American Astronomical Society 
Meeting Abstracts #216, #300.03-+ 
Miiller, D., Hansteen, V. H., & Peter, H. 2003, A&A, 411, 605 
Miiller, D., Peter, H., & Hansteen, V. H. 2004, A&A, 424, 289 
Parker, E. N. 1978, ApJ, 221, 368 
Parker, E. N. 1983, ApJ, 264, 642 
Parker, E. N. 1988, ApJ, 330, 474 
Patsourakos, S. & Klimchuk, J. A. 2006, ApJ, 647, 1452 
Peter, H., Gudiksen, B., & Nordlund, A. 2004, ApJ, 617, L85 
Peter, H., Gudiksen, B., & Nordlund, A. 2006, ApJ, 638, 1086 
Pizzolato, N., Maggio, A., Micela, G., Sciortino, S., & Ventura, P. 2003, A&A, 
397, 147 

Rappazzo, VeUi, Einaudi, & Dahlburg. 2007, ApJ, 657, L47 



S. Bingert and H. Peter: Intermittent heating in the solar corona employing a 3D MHD model 



13 



Rosner, R., Tucker, W. H., & Vaiana, G. S. 1978, ApJ, 220, 643 
Schrijver, C, Sandman, A. W., Aschwanden, M. J., & DeRosa, M. L. 2004, ApJ, 
615,512 

Solanki, S. K., Rueedi, I., & Livingston, W. 1992, A&A, 263, 339 
Solanki, S. K. & Steiner, O. 1990, A&A, 234, 519 

Spitzer, L. 1962, Physics of fully ionized gases (Interscience, New York (2nd 
edition)) 

Spitzer, L. & Harm, R. 1953, Phys. Rev., 89, 977 

Steiner, O. & Murdin, P. 2000, in Encyclopedia of Astronomy and Astrophysics, 

ed. Murdin, P. (Grove's Dictionaries) 
Ulmschneider, P., Priest, E. R., & Rosner, R., eds. 1991, Mechanisms of chro- 

mospheric and coronal heating (Springer- Verlag, Berlin) 
Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 
Warren, H. P, Kim, D. M., DeGiorgi, A. M., & Ugarte-Urra, 1. 2010, ApJ, 713, 

1095 

Withbroe, G. L. & Noyes, R. W. 1977, ARA&A, 15, 363 



